function [ DrM ] = getDrM( r, M )
%GETDRM Summary of this function goes here
%   Detailed explanation goes here

down = (r(1)^2+r(2)^2)^2;
up1 = r(1)^2-r(2)^2;
up2 = r(2)^2-r(1)^2;

s = 1;
A = M;
ds = [0, 0];
dA2 = [2*r(1)*r(2)^2, r(2)*up2, -2*r(1)^2*r(2), r(1)*up1; r(2)*up2, -2*r(1)*r(2)^2, r(1)*up1, 2*r(1)^2*r(2)] / down;
dA = reshape(dA2, 2, 2, 2);

DrM = tensorProductGradient(s, ds, A, dA);

end

